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ABSTRACT 

We have observed 99 mid-infrared-bright, massive young stellar objects and compact 
Hii regions drawn from the Red MSX source (RMS) survey in the J=3—2 transition 
of and ^^CO, using the James Clerk Maxwell Telescope. 89 targets are within 
6kpc of the Sun, covering a representative range of luminosities and core masses. 
These constitute a relatively unbiased sample of bipolar molecular outflows associated 
with massive star formation. Of these, 59, 17 and 13 sources (66, 19 and 15 percent) 
are found to have outflows, show some evidence of outflow, and have no evidence of 
outflow, respectively. The time-dependent parameters of the high-velocity molecular 
flows are calculated using a spatially variable dynamic timescale. The canonical cor¬ 
relations between the outflow parameters and source luminosity are recovered and 
shown to scale with those of low-mass sources. For coeval star formation we find the 
scaling is consistent with all the protostars in an embedded cluster providing the 
outflow force, with massive stars up to ^30 Mq generating outflows. Taken at face 
value, the results support the model of a scaled-up version of the accretion-related 
outflow-generation mechanism associated with discs and jets in low-mass objects with 
time-averaged accretion rates of Mq yr“^ onto the cores. However, we also sug¬ 

gest an alternative model, in which the molecular outflow dynamics are dominated by 
the entrained mass and are unrelated to the details of the acceleration mechanism. We 
find no evidence that outflows contribute significantly to the turbulent kinetic energy 
of the surrounding dense cores. 

Key words: starsiformation - stars:protostars - stars:abundances - stars:massive - 
stars:winds, outflows 


1 INTRODUCTION 


The formation scenario for massive protostars is compara¬ 
tively unc ertain compared with the low-mass star-formation 
paradigm dShu et al.lll987l ). During their formation, massive 
stars (>8Mq) deposit copious amounts of energy into the 
interstellar medium (ISM) and their natal cloud through 
jets, outflows and expanding ionisation fronts. These pro¬ 
cesses could act to regulate infall, accretion and local 
star formation itself. The well studied massive prot ostars 


Cep A HW2 dPatel et al.l 20 051: Curiel et al 
IRAS 20126-t4104 Jcesaroni et al' 


I 2 OOOI : ICesaronill2005l : ICesaroni et al.ii2014^ . for example, are 


200d) and 

I 999 I: ^epherd et al.l 


* drawn from the RMS survey, 

http: //rms. leeds.ac.uk/cgi-bin/public/RMS_DATABASE.cgi 
I E-mail:maud@strw.leidenuniv.nl (LTM) 


thought to be surrounded by Keplerian discs and have jets 
powering massive outflows in a direction perpendicular to 
the disc plane. The formation scenario for these particular 
sources appears to be analogous to that of low-mass proto¬ 
stars, although it is unclear whether such sources are repre¬ 
sentative of all massive protostars, or are isolated cases. 


One early signpost of star formation and a po¬ 
tential way to investigate accretion is the phenomenon 
of bipolar molecular outflows e r upting from dens e, dark 
clouds (see [Richer et al.l I 2 OO 0 I : I Arce et 'al] I 2 OO 7 I . for a 
mo re th or ough review). Early observations summarised 
in Ibadal (Il985l ) indicate the vast luminosity range of 
sonrces which drive such flows. Subsequent works have 
confirmed that cores cont aining massive protosta r s drive 
powerful outflows (e.g. She 2 herd^_Churchwell|_ 1996^ 
IZhang et al.l 1200 ll : iRidge fc Moord " 2001 : Beuther et alJ 
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2 L. T. Maud et al. 


l2002d : iKlaassen fc Wilson! l2007l : iMottram fc BruntI I 2 OI 2 II . 
The outflow phenomenon thereby provides one potential link 
between low- and high-mass star formation scenarios. 

Although detection of outflows alone is insufli- 
cient to draw conclusions on the existence of a sin¬ 
gle star-formation process, correlations between outflow 
energetics and bolometric luminosity over several or¬ 
ders of magnitude have been interpreted as evidence of 
a common outflow-driving mechan i sm that scale s with 
luminosity ll^ driguez et al.l Il 982l: iBal lv fc La d^ Il983| : 


ICabrit fc Bertoii^ 19921: Shepherd fc Churchwelll 1996ah. 

Few previous studies of outflows from massive proto¬ 
stars have used representative and unbiased samples that 
are large enough to deal with the o bservational uncertain - 
ties and intrinsic scatter in properties. iBeuther et al1(l2002cl i 
find outflows in 21 of 26 targets and see evidence of similar 
trends in a sample based on a catalogu e of massive star for¬ 
mation sites using IR AS point sources dMolinari et al.iri996l : 
ISridharan et al.ll2002l l. Such a sample, containing sources at 
a range of distances, may be subject to distance biases that 
overwhelm the significance of the correlation and suffer from 
confusion in the luminosity valu es due to the low resolu tion 
of IRAS (2-5 arcmin at 100 /rm). iRidge fc Moord (120011 1 ex¬ 
amined a sample of 11 same-distance (2-kpc) sources and 
found only tentative evidence for a scaling of outflow en¬ 
ergetics with bolometric source luminosity in the massive 
prot ostar regime, althoug h their sample was small. 

iDuarte-Cabral et al.l ll2013ll investigated 9 outflows 
from mid-IR-dark, high-mass protostellar analogues of 
Class-0 low-mass sourcesQat ~l-2 arcsecond resolution and 
found that the outflow pr operties scale over '^2. 5 orde rs of 
magnitude in luminosity. ISan Jose-Garcia et al.l (l2013ll ex¬ 
amined the CO properties of a sample of low, intermediate 
and high-mass sources, they found a general trend of in¬ 
creasing outflow velocities with bolometric luminosity, but 
with significant intrinsic scatter. Still, larger samples are re¬ 
quired to investigate outflow dynamics and kinematics from 
representative selection of massive protostel lar regions. 

The Red MSX Source (RMS) survey llLumsden et al.l 
I200II2OI3I) specifically identifies massive young stellar ob¬ 
jects (MYSOs) and Hll regions drawn f rom the MSX mid - 
infrared survey of the Galactic plane (lEgan et al.l l2003ll , 
which has higher angular resolution (18 arcseconds) than 
earlier IR AS point-source catalogues of massive star fo rma- 
tion sites llMolinari et al.|[l996l : l^idharan et al.l 1^021 ! and 
therefore is less affected by confusion in the Galactic Plane. 

This is the second of two papers investigating the 
same sample of 99 MYSOs and very compact Hll regions 
drawn from the RMS sample. It is the first large sample 
of molecular outflows from high-mass and high-luminosity 
YSOs that is not dominated by malmquist-type biases. 
In many previous attempts to loo k at high-mass YSOs 
(e.g. IShepherd fc Ghurchwelll Il996al : IBeuther et al.l l2002d . 
and the high-mass, >1 0^ Lq component of the definitive 
iGabrit fc Bertouall992l study), accounting for the distance 
dependencies and possible biases reduces the significance of 
the correlations. The conclusion of a single driving mecha- 


^ youngest observed YSOs that are deeply embedded 
(Metm) >Matar) an d h ave l arge sub-mm excesses, see 
iBontemps et al.l lll996l l and iLadal lll987h 


nism in both low- and high-mass outflow sources therefore 
requires more careful examination^ _ 

In the first of the two papers, iMaud et al.l (l2015ll (here¬ 
after Paper I), we identified and analysed the properties of 
the dense cores in these massive star-forming regions. The 
luminosities and masses of the cores were shown to be un¬ 
affected by significant distance biases out to 6kpc, and are 
representative of the range of luminosities in the whole RMS 
MYSO/Hll-region sample subject to a similar distance cut 
(~250 objects). Note we use the terminology ‘cores’ through¬ 
out this paper in analogy with other studies of similar high- 
mass-star formation regions with sim ilar resolutions (see. 
IBeuther et al.l l2002bl : I Hill et ahl 120051 ') . We probe a range 
of radii for these dense (>10^cm“®) multiple star-forming 
‘cores’ (-^ 0.1 — 1.0pc), with a nominal value of 0.35pc 
(Paper I). As noted in paper I, we understand tha t there 
is further sub-structure (e.g. iBontemps et al.ll2010lj which 
may pertain to single s tar-forming cores much smaller than 
~ 0.1 pc in radius (e.g. iHennemann et al.|[2009ll and all our 
cores are interpreted as containing young stellar clusters, 
in which many lower-mass protostars are likely to be co¬ 
evolving with the high-mass objects. Furthermore, in terms 
of a mass-luminosity plot, these regions are indistinguishable 
from the larger sample; i.e., the MYSOs and Hll regions are 
at a similar evolutionary stage (also s ee Davies et al.ll201ll : 
iMottram et al.l[2011al : lUrauhart et al.ll2014bh . 

Here the JGMT HARP and ^^GO (J=3-2) ob¬ 

servations are specifically used to identify and examine out¬ 
flows emanating from these regions. Its higher critical den¬ 
sity and upper-level energy make the (J=3—2) transition an 
ideal tracer for investigating outflows when compared with 
the (J=l—0) transition, as the entrained outflow emission 
is warmer than the core and because line-of-sight ambient 
emission is less prevalent, therefore less likely to contam¬ 
inate outflow wings. The distance-independent nature and 
large sample size allows a definitive investigation of outflows 
from massive protostars. 

Section [ 2 ] summarises the sample and the observations 
undertaken. Section |3] describes the outflow identification, 
parameter determination and outflow parameter calcula¬ 
tions. Section |4] presents the general results obtained and 
Section [5] discusses the observed trends, comparisons with 
lower-mass outflows, the outflow-driving mechanism and the 
effect of the outflows on the natal cores. The main results 
are summarised in Section!^ 


2 SAMPLE AND OBSERVATIONS 

Paper I provides a more comprehensive overview of the sam¬ 
ple and selection criteria. To summarise here, the 99 sources 
are all MYSOs and compac t Hll regions from the RMS 
sample (iLumsden et al.ll2013l l. The distance limited statis¬ 
tical sample corresponds to 89 of these sources at heliocen¬ 
tric distances less than 6kpc. Except where indicated, the 
luminosities of all sources w ere calculated using the most 
up-to-date source distances llUrauhart et al.112012| . l2014a^ 
and multi-wavelength spectr al energy distribution (SED) fits 
from IMottram et ahl (l2011bf l. 

The observations were undertaken using the James 
Glerk Maxwell Telescope (JGMT) in 2007 and 2008 as part 
of projects 07AU08, 07BU16, 08AU19 and 08BU18. The full- 
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Massive molecular outflows 3 


width half-maximum (FWHM) beam size at ~345 GHz for 
the (J=3—2) transition is arcsec. The Hetero¬ 

dyne Ar ray Receiver Progr am (HARP) 16-pixel SSB SIS 
receiver dBuckle et al.llioogl l was used with the ACSIS cor¬ 
relator (Auto-Correlation Spectral Imaging System) back¬ 
end configured with an operational bandwidth of 1000 MHz 
for the transition. For the data, taken simultane¬ 

ously with C^®0 (see Paper I), the bandwidth was 250 MHz. 
The effective velocity resolutions were -^0.4 (^^CO) and 
~0.06kms“^ (^®CO and C^®0). In calculations of the op¬ 
tical depth (see Section o the data is re-sampled to 

match the velocity resolution. Typical spectral noise 

levels (dTmb) for both the and data range be¬ 

tween -^0.4 and 0.6 K in a ~0.4kms“^ bin. 

Maps of the sources were taken in raster-scan mode with 
continuous (“on-the-fly”) sampling and position switching 
to observe a ‘clean’ reference position at the end of each 
scan row. Pointing was checked against a known, bright 
molecular source prior to each science observation and is 
likely to be within ~5arcsec, typical of the JCMT. Reduc¬ 
tion was undertaken with a custom pipeline utilising the 
KAPPA, SMURF, GAIA and SPLAT packages which are part of 
the STARLINK software maintained by the Joint Astronomy 
Centre (JAC)0. Linear baselines were fitted to the source 
spectra over emission-free channels and subtracted from the 
data cubes. Bad baselines on each of the working receivers 
were flagged if they showed evidence of sinusoidal fluctua¬ 
tions. The final cubes used in the analysis were re-gridded 
to a 7-arcsec spatial pixel scale. The data, originally on the 
corre cted antenna temperature scale (TX; iKutner fc Ulichl 
Il98ll l were converted to main-beam brightness temperature 
Tmb = TX/rjmb, where T^mb = 0.6 6 as measured by J AC dur¬ 
ing the commissioning of HARP (iBuckle et al.11200^ 1 and via 
continual planet observations. 

Table [T] lists the sources and their properties as ex¬ 
tracted from the RMS online data base0. In addition, other 
common source association names are presented along with 
the IRAS designation and offset if available. 


3 OUTFLOW IDENTIFICATION AND 

PARAMETER DETERMINATION 

Outflows are typically identified by characteristic high- 
velocity wing components in molecular rotational emission¬ 
line spectra. Outflows with axes close to the plane of 
the sky will not exhibit such broad wings but may have 
a linear structure in an integrated emission map encom¬ 
passing near-ULSR velocities. However, even a known low- 
inclination, high-mass outflow source in onr sample, GGD- 
27 (HH80-81, GOIO.8411—02.5919), has reasonably broad 
wings with -^15kms~ ^ full-width-zero-intensity (FWZI) 
dYamashita et aT1ll989l l. 

The data were investigated interactively (using the GAIA 
package), by examining slices of the data cube while simul¬ 
taneously extracting spectra within a 3-pixel diameter cir¬ 
cular aperture. To highlight the outflow emission, a 


^ http://starlink.eao.hawaii.edu/starlink 

® http://rms.leeds.ac.uk/cgi-bin/public/RMS_DATABASE.cgi 


cut is applied to the cubes (where is the standard de¬ 
viation of the spectral noise measured in emission-free ve¬ 
locity channels extracted at every pixel). Figure [T] illustrates 
the identification process for the sonrce G078.1224-1-03.6320 
(IRAS 20126+4104), a well stndied MYSO with a bipolar 
ontflow. Slices of the cube at blue- and red-shifted velocities 
with respect to the source ulsr are shown and inset are the 
spectra extracted from the cube at spatial positions corre¬ 
sponding to the white circle. It Is clear that the blue-shifted 
outflow component lies to the north while the red- shifted 
material lies to the south (c.f. ISheoherd et al.l|2000fl . Note 
the more chaotic spatial structure closer to the ulsr veloc¬ 
ity (—3.6kms~^). This investigation was undertaken for all 
targets in order to establish the presence or lack of extended 
velocity material over the core extent found in Paper I as a 
proxy for an outflow. 


3.1 Temperature Calculation and Optical Depth 
Correction 

The properties of the outflows, mass, momentum and energy 
are determined, essentially, from the moments: 

III N{x,y,v)v" dxdydv [n = 0,1, 2] (1) 

all of which hinge on an accurate determination of the col¬ 
umn density N{x,y,v) (where x,y are spatial coordinate 
and V the velocity coordinate in a data cube). To estimate 
the column density elements N{x,y,v), the ^^GO (3—2) 
emission is assumed to be in local thermodynamic equi¬ 
librium (LTE). Paper I notes that this is likely to be the 
case for these cores. A single excitation temperature Tex and 
position- and velocity-dependent optical depth values r are 
calculated for each source. As the '^^GO emission is shown 
to be optically thick (at the core ulsr) in all but four of 
these sources (Paper I), it is expected that the ^^GO is also 
optically thi ck. This Is known to be the case in other o utflow 
sources (e.g. iGabrlt fc Bertou3ll990l : IChol et ^Il993ll . 

Thus a single Tex for each source is calculated from: 


16.59 K 

ln[l + 16.59 K/(r^b,i 2 + 0.038)] ’ 


( 2 ) 


where hiyC^CO)/k = 16.59 K, with uC^CO) = 345.80 GHz 
and rmb,i 2 is the peak main-beam brightness temperature of 
the emission at the location of the source within one 

beam and at the source ulsr- Note, in some cases where self 
absorption occurs the velocity of the peak rinb,i 2 is slightly 
shifted from the source ulsr- 

Often, a constant factor is used to correct outflow 
masses for optical depth (e.g. iBeuther et al.l[2002cl l. How¬ 
ever, since we have both and data cnbes, we 

can calculate an optical depth for each voxel and correct ev¬ 
ery spatial {x, y) and velocity (u) element by its own optical 
depth, assuming Tex is the same for both isotopologues, in 
which case: 


r^b,i2 ^ 1 - ^ 1 - 

Tnib.is 1 - e-^13 1 - ’ 

where R is the abundance ratio of [^^CO]/[^^CO] = 7.5 x 
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Table 1. Source parameters for all objects in the sample, taken from the RMS survey online archive. Only a small portion of the data 
is provided here, the full table is available in the electronic supplementary information and via the RMS database directly. 


MSX Source Name 

RA. 

(J2000) 

DEG. 

(J2000) 

Type 

vi.sn 
(km s“^) 

Distance 

(kpc) 

Luminosity 

(L©) 

IRAS source 
(offset) 


Other 

Associations 

GOlO.8411-02.5919 

18:19:12 

-20:47:30 

YSO 

11.4 

1.9 

24000 

18162-2048 

(4») 

GGD27 

G012.0260-00.0317 

18:12:01 

-18:31:55 

YSO 

110.6 

11.1 

32000 

18090-1832 

(3") 


G012.9090-00.2607 

18:14:39 

-17:52:02 

YSO 

35.8 

2.4 

32000 

18117-1753 

(11") 

W33A 

G013.6562-00.5997 

18:17:24 

-17:22:14 

YSO 

48.0 

4.1 

14000 

18144-1723 

(2") 


G017.6380-I-00.1566 

18:22:26 

-13:30:12 

YSO 

22.5 

2.2 

100000 

18196-1331 

(11") 


G018.3412-I-01.7681 

18:17:58 

-12:07:24 

YSO 

32.8 

2.9 

22000 

18151-1208 

(16") 


G020.7438-00.0952 

18:29:17 

-10:52:21 

Hii 

59.5 

11.8 

32000 



GRS G020.79-00.06 

G020.7491-00.0898 

18:29:16 

-10:52:01 

Hii 

59.5 

11.8 

37000 



GRS G020.79-00.06 

G020.7617-00.0638 

18:29:12 

-10:50:34 

YSO/Hii 

57.8 

11.8 

62000 



GRS G020.79-00.06 

G023.3891-I-00.1851 

18:33:14 

-08:23:57 

YSO 

75.4 

4.5 

24000 

18305-0826 

(6") 

GRS G023.64-I-00.14 


Dgc + 7.6, following IWilson fc R.oodI (|l994h . and Dgc is 
Galactocentric distance in kpc. Equation [3] can be used ef¬ 
fectively where the emission is optically thin, and al¬ 

though -fia^core ~ 3.5 (Paper I), ria <1 in the line wings 
aw ay from the so urce wlsr. 

iMaudl (l2013h found that the mean of all the individ¬ 
ual voxel {x, y, v) optical depths, ri 2 , is consistent with that 
derived using !](/ T’mb,i2 dw) and !](/Embps dw) directly 
in place of rmb,i 2 and Tmb.ia in Equation |3] However, later 
calculations of momentum and energy, for example, can be¬ 
come grossly overestimated (by orders of magnitude) using 
such a single, averaged optical-depth correction, because the 
higher velocity, optically thin emission is over-corrected. 


3.2 Outflow Masses 

Once the optical depth and excitation temperature are es¬ 
tablished for each voxel, the column density maps for each 
outflow lobe, N{x,y) pixels, are calculated via: 


iV(i2cOx,y) = 4.78 X 10^^ e^P(16.74/rex) (Tex + 0.93) 


/ 


Tmb 


exp(-16.59/Tex) 
_ n2 _ 

[1 - exp(-ri 2 )] 


dv , (4) 


where Af(C^^O) is in cm“^ and dv is in this case is the 
integration over the velocity extent in kms“^; hence, the 
mass elements at each pixel coordinate, x, y (where n = 0 
in Equation [TJ are. 




iV('^COx,y) 


■ H2 
12CO 


Ug D'^ 


( 5 ) 


where Q is the solid angle of a pixel element (x,y), D is 
the distance to the source, (H 2 /^^CO) is the H 2 to ^^CO 
abundance ratio = 10"^ and /ig = 1.36 is the total gas mass 
relative to H 2 . A more detailed derivation of the column 
density and mass are given in Appendix lAl 

Calculation of the outflow mass, which then leads to mo¬ 
mentum and energy (where n = 1,2 in Equation[T]), however, 
first relies on the correct identification of emission associated 
with the high-velocity gas in order to reliably separate it 
from that of the core. The method used here is analogous to 
that outlined in Paper I, a combined process of integration 


over velocity (i.e. creation of zeroth-order moment maps for 
each outflow lobe) followed by an aperture summation over 
the area of outflow emission. The blue and red-shifted out¬ 
flow velocity extent (maximal velocities) are established via 
a direct, manual, investigation of the raw data cubes, chan¬ 
nel by channel, until emission spatially associated with the 
outflow drops below 3CTT„b within a beam area. Automating 
this process for a large sample is difficult as there is typically 
diffuse emission at ‘outflow’ velocities elsewhere in the maps 
that is not associated with the outflow itself. 

iRidge &: Moord (1200 il l note that careful removal of the 
core component is required to determine outflow properties 
accurately. The low-velocity limits of each outflow are es¬ 
tablished with reference to the velocity extent of the core 
component in Paper I. We use the largest of; (1) the es¬ 
tablished velocity range (Paper I, Table 3); or (2) the full- 
width-tenth of maximum (FWTM = 1.8 x C“0 FWHM) 
as this ensures that most of the core emission is excluded 
from that associated with the outflow lobes. Generally, case 
1 ~ case 2 wherever the G*^®© emission is well above 
However, this is not the case for all the cores due to varying 
noise levels; hence the rationale of using the maximum of 
case 1 or 2. Furthermore, the difference between case 1 and 
case 2 is not more than 1.0 kms“^, which is accounted for in 
the uncertainty analysis. At most, we find the mass of the 
outflow could be underestimated by 50 percent. However, 
the influence on the momentum and energy parameters is 
smaller as the ‘missed’ emission is at very low velocities and 
contributes much less to values containing multiples of the 
velocity. The velocity ranges established by this process are 
listed in Tabled including the spectral noise, and the 

maximum velocities shifted by the source ulsr- 

After integration over the selected velocity range, a 
polygon aperture is used to sum the spatial extent of the 
outflow emission in the integrated lobe maps. A problem 
with this method is that true plane-of-the-sky outflows (i.e. 
with no wing components) will be excluded because their 
velocity range is wholly within that attributed to the core. 
However, we do not identify any plane-of-sky outflows in the 
cube investigation stages, where we search the data cubes 
directly for evidence of elongated features close to the core 
ulsr- 

The polygon aperture used to define the spatial extent 


© 2014 RAS, MNRAS 000, [T]-?? 






















Massive molecular outflows 5 




Figure 1. Each sub-plot shows a velocity slice (indicated in top left of each sub-plot) from the data cube for G078.1224+03.6320. Inset 
is the summed spectrum extracted within the white circle. It is clear when moving from blue to red-shifted velocities that there is an 
outflow with a clear spatial offset about the central source location (white cross). 


of an outflow is adapted from the simple 3 ctmap contour 
method established in Paper I. The ctmap is calculated from 
the average spectral noise, and is equivalent to a spatial av¬ 
erage of the physical map noise (see Appendix B), but we 
And that only in a handful of cores, with little to no ambient 
emission in the regions surrounding the outflow lobes, 
does 3 Ctmap correctly delineate the outflow lobe region (e.g. 
as in G078.1224+03.6320). In more complex regions, with 
ambient emission away from the main outflow, we can use 
the median value of the map to act as a background level 
(note the median map values are zero, for sources with no 
ambient emission, such as G078.1224+03.6320), and define 


the aperture delineating the outflow lobes to follow the con¬ 
tour of 3 Ctmap above this background (see Appendix B for 
details). 

However, we note that, in some sources the outflow ar¬ 
eas are not well defined by this method and the apertures 
do have to be manually adjusted to divide the diffuse emis¬ 
sion from that of the outflow. These sources are flagged as 
‘selected’ (S) apertures in Table (2] Appendix C presents the 
outflow lobe maps for each source and the apertures used 
for area summation, while Figure [5] shows the example of 
G078.1224+03.6320. Note that apertures are set prior to 
optical depth correction as they are based on the observa- 
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Table 2. Outflow-detection parameters for all objects in the sample. Y, M, N in the flow column indicates whether the source has an 
outflow, shows some evidence of an outflow (maybe), or no outflow at all. Y, S, N in the Aperture Flag column indicates a good aperture, 
a manually ‘selected’ aperture or no aperture. Sources without outflows have no aperture, but also some M sources have no aperture 
due to complex diffuse emission. Av is the raw velocity extent of the lobes with respect to the observed velocities, whereas Umax are the 
maximum velocity offsets from the ulsr- The full version of this table is available in the electronic supplementary information. 


MSX Source Name 

Flow 

Flag 

Spec. 

Noise 

(K) 

A'Wblue 

(kms“l) 

(k 

ms-l) 

Blue Map 
Med. Noise 
(Kkms“l) 

Blue Map 

(cr) 

(Kkms“l) 

Red Map 
Med. Noise 
(Kkms“l) 

Rod Map b I'max.r 

(ct) 

(Kkrns”'^) (kms“'^) (kms“l) 

Aper. 

Flag 

GOIO.8411-02.5919 

Y 

0.4 

( -6.1 

, 10.0 ) 

( 14.2 

28.2 ) 

0.1 

1.2 

3.4 

1.1 

18.2 

16.0 

Y 

G012.0260-00.0317 

Y 

0.9 

( 93.9 ,108.1 ) 

(113.8 ,123.1 ) 

6.7 

2.2 

1.0 

1.8 

17.0 

12.2 

Y 

G012.9090-00.2607 

Y 

0.5 

( 7.0 

, 32.3 ) 

( 40.4 
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Figure 2. Left: Map of integrated intensity in the blue-shifted outflow lobe of G078.1224-1-03.6320, prior to optical depth correction. 
The integration velocity range is —40.0 to —5.7kms“^. The blue dashed line indicates the aperture at the contour equivalent to 3 (Tmap 
above the background (defined by the median map value). The black contours indicate the 30, 50, 70 and 90 percent levels to highlight 
offset of the emission from the source location indicated by the white cross. Right: As left but for the red-shifted outflow lobe in the 
range is —1.0 to 40kms~^. The source ulsr is —3.6kms“^. The source name, distance and outflow flag are indicated at the top of the 
plot, Y, in this case represents a good outflow source with well defined apertures. 


tional noise levels. All the appendices are available in the 
additional online material. 


3.3 Outflow Parameters 

Empirical tests conducted as part of the analysis investi¬ 
gated a number of methods u sed in the lit erature to establish 
the outflow parameters (see iMau 3 120131. for more details). 
The most accurate in calculating the momentum, P, and 
energy, E, are: 

P = Y (6) 

x,y,i 


P = i Y (7) 

x,y,i 

where M is the mass, i represents each velocity bin, at ve¬ 
locity Vi with respect to the source wlsr summed over all 
velocities and spatial pixels for both outflow lobes (as Equa¬ 
tion [T] with n = 1 and 2). These are the most accurate in 
recovering input parame ters when tested on outflow models 
dCabrit fc Bertoud Il99flh . This method of full cube analy¬ 
sis (x, y, i) avoids the overestimate of parameters that oc¬ 
curs when the total outfl ow mass is simply multiplied by the 
maxi mum velocity (e.g. ISeuther et al.ll2002cl : iLebron et al.l 
I 2 OO 6 II . a method that noti onally places all the m ass at the 
maximal outflow velocity dMargulis fc LadallQSSll . 
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In order to establish the mass flow rate, Mout, the force 
(momentum supply rate), P = Fm, and the mechanical lu¬ 
minosity (outflow power), E = Lm, a dynamical timescale 
fdyn = R/v must be estimated, where R is the distance of 
the outflow from the source. The dynamic timescale is com¬ 
monly taken to be an indicator of the age of the outflow 
and, as s uch, could also be rep resentative of the protostar’s 
age isee. iBeuther et aDl2002a i (under the assumption of a 
constant outflow over time that begins as soon as the proto¬ 
star is form ed). Further interpre tation is discussed below in 
Section 15.21 iLada fc Fichl (Il996l i introduced a means to cal¬ 
culate tdyn at all spatial positions using Rx,y/{'Ux,y), where 
{Vx,y) is the intensity-weighted-mean outflow lobe velocity, 
representative of the bulk motion of material, and calculated 
from {Px,y/Mx,y)■ To utilise our outflow lobe maps we adopt 
this method. From the position-dependent tdyn, the dynamic 
parameters are calculated via Maut(x,y) = M^{x,y)l'tdyn(x,y), 
R(x,y) Rm(x,y) R{x,y)/^dyii{x,y) aud E(^x,y) ^ 111 ( 3 :,y) 

E(x,y)/tdyn(x,y)- Ih Figure [3] it is clear that the largest con¬ 
tributor to M, P and E are the spatially offset, higher- 
velocity outflow components. The spatially diffuse, chaotic, 
low-velocity emission contributes very little (~ a few per¬ 
cent) to the total values. 

From all outflow parameters the mass is the only one 
independent of velocity. The velocity dependent variables 
(and those involving outflow length), in principle, should be 
corrected for the source inclination angle as we are only sen¬ 
sitive to the line-of-sight velocity component (or that mea¬ 
sured on the sky, in terms of lengths). In light of this: the mo¬ 
mentum and energy will be -slower limits (both oc 1/cos" 9, 
where n = 1 and 2 respectively, and where 6 is the incli¬ 
nation angle with respect to the line of sight); the dy¬ 
namical timescale will be over- or under-estimated if 9 is 
> 45° or < 45° respectively (oc 1/tan 9)-, Mout follows the in¬ 
verse relationship of tdyn; and Pm and Lm will be over- or 
under-estimated where 9 is less-than or greater-than -^38° 
(oc sin^/cos" 9 where n = 2 and 3). Inclination angles are not 
available for these sources and cannot easily be established. 
Ivan der Marel et al.l (l2013l l apply inclination angle correc¬ 
tions for their sample of low-mass outflows, using the no¬ 
tion of pole-on, plane-of-sky, and somewhere between; often 
however, a constant corr ection factor is used fo r a nominal 
inclination of ~57.3° fe.g iBontemps et al.lll99^ f. We do not 
apply any factors unless specifically noted (see Sections 15.II 
and 15.51) . 


4 RESULTS 

From the total sample of 99 sources 65 are associated with 
molecular outflows, all of which exhibit line-wings in their 
spectra. 14 show no evidence of outflow identifiers 
such as high-velocity wings or plane-of-the-sky linear struc¬ 
tures - their and spectra are purely Gaussian in 

profile. 20 sources in complex regions (spatially and spec¬ 
trally) show some high-velocity components but also multi¬ 
ple spectral peaks which confuses the outflow identification. 
In the distance-limited subset (89 sources) there are 59 def¬ 
inite outflows, 17 with outflow like properties and 13 with 
no outflow evidence. Not all of the sources with outflows 
show a clear spatial offset of the blue and red-shifted veloc¬ 
ity components like G078.1224+03.6320. Given the resolu¬ 


tion of the single-dish observations, it is not expected that 
all spatially offset velocity components would be resolved. 
There is evidence, however, for a distribution of source in¬ 
clinations, as some of the most distant sources, expected to 
be the least resolved, do exhibit clear spatially offset blue- 
and red-shifted outflow lobes (e.g. G053.5343—00.7943 at 
5 kpc). Only two sources appear to have clear outflows close 
to the plane of the sky associated with the cores (the afore¬ 
mentioned GOIO.8411-02.5919 and G109.8715+02.1156, al¬ 
though these cores drive multiple distinguishable outflows, 
see Figure mi. This can also be confirmed as the sources 
have been previously well studied at radio wavelengths and 
have very line ar radio jets (GGD 27. iMarti et al.l[l993l and 
Cep A HW 2 . ICuriel et al.ll 200 ^'l . G203.3166+02.0564 also 
has a close to plane-of-sky outflow forming a linear struc¬ 
ture across the map. This outflow is offset from the source 
location to the East and the velocity range is optimised only 
for the outflow ass ociated with our source position (also see 
iMaurv et al.l [ioo^ . Cunningham et al. 2015, submitted to 
MNRAS). However, as previously noted we do not identify 
any isolated, true plane-of-the-sky outflows. It is inherently 
difficult in these sources to spatially separate the complex 
core emission from potential plane-of-the-sky outflows and 
so it is possible that some may be missed as they are simply 
confused by the core emission. 

Table [3] lists the masses, momenta and energies while 
Table 0] lists the dynamical-time-dependent parameters, 
mass flow rate, force and power (luminosity). The reported 
values of Mout, P and E are those summed within the de¬ 
fined apertures (blue and red dashed contours in Figures [2] 
and 0 , although the tdyn values in Table 2] are those cal¬ 
culated directly from i?max/umax, as one cannot establish a 
meaningful single value from a spatially variable map (note 
the values obtained for Mout, P and E using this single tdyn 
are closely comparable with those using presented in Table 
fusing the spatially variable t^yu(x,y), see Section [5. 51) . 

Figure [S] shows the dependence of the calculated out¬ 
flow parameters M, P, E, M, P and E on the bolometric 
luminosity of the associated MYSO or Hll region in the RMS 
catalogue (shown in Table [T|). These plots also include the 
few sources where D> 6 kpc, although they are not used in 
the statistics or analysis. Table [5] provides a list of all Spear¬ 
man rank correlation statistics for relationships in Figure 
(5] using only the sources in the distance-limited subsample. 
Table [ 6 ] gives the parameters of the corresponding linear fits 
to the data, again using only the distance-limited sources. 


5 DISCUSSION 

5.1 Mass, Momentum and Energy 

The outflow mass is the only velocity-independent variable 
and so the only one not affected by source inclination an¬ 
gle. The top-left panel of Figure [5] indicates that the cores 
harbouring the most luminous protostars are those with the 
most massive outflows. Both the blue- and red-shifted out¬ 
flow lobe masses are plotted to show that both follow the 
same trend, and to investigate any possible systematic trend 
or differences between the lobes. We continue to plot both 
lobe properties thought the analysis. Even though asymmet¬ 
ric or even single-lobed outflows are not uncommon, in most 
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Figure 3. Top—Bottom: The spatial distributions of the calculated dynamical timescale, mass flow rate, force and outflow power 
(luminosity). The left and right panels show the blue- and red-shifted outflow lobes, respectively. With reference to Figure [2] the 
positionally offset, higher-velocity emission corresponds to lower dynamical timescales. The mass flow rate, force and power are strongest 
in regions of low tdyn(a:,y )5 expected. Lower velocity diffuse material in the case of G078.1224-1-03.6320 contributes little to the flow, 
force and power as these have the largest tdyn(a:,i()- Note the logarithmic scaling for flow and force to show the weak emission. 


© 2014 RAS, MNRAS 000, [T]-?? 























































Massive molecular outflows 9 


G109.8715+02.1156 0.7 kpc 


o 

0) 

w 

o 

O 

LU 

Q 

0) 

w 


100 

' 1 'j..'.-' 1 .j 1 1 1 1 1 ' ' ' ' 1 ' ' ' ' 1 ' ' ' i 

\ :> 

r. '•! 

S 1 


' 1 ' ' ' ' 1 ' ' ' ' 1 ' ' ' ' 1 ' ' ' ' 1 ' ' ' I- 

50 

;■ ‘i. 



0 

Is' :' 



-50 



■ " \ '^-^1 : 

100 



/• ... 

\ •, 

- .f 
-I 

. 1 . . . . 1 . . . All . . .^'.'^1 .... 1 


100 50 0 -50 -100 

offset RA (arcsec) 

100 50 0 -50 -100 

offset RA (arcsec) 


Figure 4. As Figure[2]but for the source G109.8715+02.1156. The velocity ranges are listed in Table[2] This source contains at least two 
outflows, one close to the plane of the sky (extension in a roughly North-South direction), and another, with broader outflow wings, in 
an East-West direction. Given the spatial resolution, separation of these components is subjective and hence the total ‘cluster’ outflow 
parameters are calculated (see Section l5.5ll . 


cases where both blue- and red-shifted lobes are present, the 
masses are consistent within a factor of '^2—3. This is con¬ 
sistent with what one would expect if outflows are entrained 
material in an isotropic core, for example. A few sources 
have very asymmetric profiles and, in two particular cases 
(G023.7097+00.1701 and G050.2213-00.6063), there is an 
absence of high-velocity emission in the blue and red lobes, 
respectively. There is no clear preference for more massive 
blue- or red-shifted outflows in the sample. As discussed in 
Paper I, these cores are clusters of protostars, each of which 
could be powering outflows that combine and interact such 
that we observe a preferentially stronger blue- or red-shifted 
flow. Furthermore, dependent on the core geometry or den¬ 
sity distributions the material constituting to the outflow 
itself may vary considerably (also see Section [5.211 . 

The next two panels (top right and left middle) of Fig¬ 
ure [3 show the same clear power law trends of increasing 
momentum and energy in more luminous protostars seen in 
many previous studies. The simplest interpretation is that 
the jet or wind from the most luminous protostar in each 
core is able to entrain more of the available core material and 
thus drive the most powerful and most energetic outflow in 
the region. Assuming a young stellar cluster is present in 
these sources, and that its luminosity is dominated by the 
most massive cluster member, then these relationships sup¬ 
port the idea of an outflow-driving mechanism that scales up 
to protostars as massive as ~50 Mq (~5x 10® Lq), our most 
massive source (including those more distant than 6 kpc). 

There is no apparent difference in mass, momentum or 
energy derived for sources classed in the RMS survey as 
YSO or Hll regions, even thou gh these are often thoug ht of 
as different evolutionary states. Urauhart et al.l ll2014bl) find 
no difference in the Lboi versus Mcore distribution for these 
two classifications. We also noted in Paper I that the core 
properties of these sources are indistinguishable and so they 


are likely to be at roughly the same evolutionary stage, with 
similar outflow properties. 

Some of the additional scatter in the momentum and 
energy values in Figure [5] is likely caused by a distribu¬ 
tion of outflow inclination angles. The spatial resolution 
of these single-dish observations does not allow inclina¬ 
tion angles t o be established . The mean inclination angle 
of ~57.3° in iBontemps et al.l lll996l) would results in only 
a constant scaling of >^1.85 (l/cos0) for momentum and 
~3.43 (l/cos^0) for energy, and therefore does not change 
the relative relationships seen in Figure[5]i.e. in log-log plots. 
ICabrit fc BertoutI lll990l ) however, discuss how outflows with 
large inclination angles, >70°, could have energies underes¬ 
timated by an order of magnitude. Furthermore, as already 
distinguished in two of our cores, there are likely multiple 
outflows driven by the sources within the cores, and hence 
a single inclination angle may not be fully representative. 

5.2 Dynamical Timescales 

The use of tdyn assumes that the local gas velocity is equal 
to the velocity of the shock wave driven through the molec¬ 
ular gas by the underlying wind or jet, which is the case 
for an isothermal shock. Sin ce the shock velocit y is likely 
to be position-dependent, the lLada fc FictJ (Il996l) position- 
variable tdyn method is more physical but does not easily al- 
low a single chara c terist ic value to be established. However, 
iDownes fc CabritI (120071 ) report that calculating i?max/t>max 
is t he best ‘classi c al’ wa y of obtaining a single tdyn value 
(as iBeuther et ahl (l2002d ). for example). Figure [S] therefore 
presents dynamical timescales, calculated via the simplistic 
Rmax/umax for the distance-limited subsample (D<6kpc), 
against luminosity, showing no significant correlation (Ta¬ 
ble [5]). 

We find that the dynamical timescale dependent pa- 
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Figure 5. Top—Left to Bottom—Right, plots of outflow mass, momentum, energy, mass flow rate, force and power against the source 
luminosity. The blue and red symbols indicate the blue- and red-shifted outflow lobe values and are joined by a dashed line for each 
source. Open and fllled symbols are MYSOs and Hii regions, respectively, while circles and squares are for sources with ‘Y’ and ‘S’ 
aperture flags. The sources that have uncertain outflow evidence and flagged as ‘M’ are indicated by the x symbols. The source symbols 
that are circled are those with D>6kpc. The correlation of all parameters with source luminosity is clear. Errors for mass, momentum 
and energy those calculated by increasing and decreasing the integration velocity range by one bin (~0.4kms~^) at the upper and 
lower limit, while those indicated for flow, force and power are 50 percent error bars. These uncertainties are direct from the change in 
integrated velocity range and do not account for uncertainties in source distanc e, or more import a ntly t he outflow inclination (except 
mass), and should be considered as minimal uncertainties. As noted in the text. ICabrit & Bertoud l|l990h indicate that uncertainties in 
the energy can be up t o an order of magnitude for large inclination angles (> 70°). The error shown for source luminosity is 30 percent 
i|Mottram et al]l2Qllblb 
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Table 3. Mass, momentum and energy calculated for all sources (including those where D>6kpc) with outflows or with evidence of 
outflows where apertures could be defined. The ‘b’, ‘r’ and ‘total’ subscripts indicate the blue-shifted lobe, red-shifted lobe and total 
values of each parameter. The mass is given in units of Mq, momentum in Mq kms“^ and energy in 10“^^ erg. Uncertainties calculated 
from adoption a different integrated velocity range (di~0.4km~^ at the upper and the lower velocity boundary) are 26 and 23 

percent for mass, momentum and energy respectively. Note, error s in source distance, and more importantly outflow inclination angle 
could have a much more significant effect. ICabrit &: BertoutI lll990h suggest uncertainties up to an order of magnitude for energy values 
if the outflows have large inclination angles (> 70°). The full table is available online. 


MSX Source Name 

Mb 

Mr 

^total 

Pb 

Pr 

Pfotal 

Eb 

Er 

Efotal 

GOIO.8411-02.5919 

21.3 

33.4 

54.7 

62.9 

97.1 

160.0 

225.2 

342.1 

567.4 

G012.0260-00.0317 

81.0 

56.4 

137.4 

440.8 

251.0 

691.8 

3201.9 

1306.8 

4508.8 

G012.9090-00.2607 

18.5 

60.9 

79.3 

121.1 

580.2 

701.3 

1068.8 

6948.0 

8016.9 

G013.6562-00.5997 

35.8 

4.3 

40.1 

174.9 

32.0 

206.9 

1005.1 

301.6 

1306.6 

G017.6380+00.1566 

56.6 

40.7 

97.3 

235.6 

224.7 

460.3 

1209.6 

1409.0 

2618.6 

G018.3412+01.7681 

48.3 

18.5 

66.8 

92.4 

81.8 

174.2 

197.2 

430.8 

628.0 

G020.7438-00.0952 

372.3 

94.1 

466.4 

1837.5 

361.8 

2199.3 

10296.9 

1490.0 

11786.9 

G020.7491-00.0898 

264.7 

94.1 

358.8 

1429.9 

388.2 

1818.1 

8438.5 

1699.7 

10138.2 

G020.7617-00.0638 

64.7 

115.6 

180.2 

284.7 

634.2 

919.0 

1359.2 

4884.6 

6243.8 

G023.3891+00.1851 

28.4 

44.7 

73.0 

108.0 

198.7 

306.8 

612.8 

1276.9 

1889.7 


rameters M, P and E, obtained using the single idyn = 
-Rmax/umax, via M ~ M/tdyn, otc., are comparable with 
those found using the positio n-variable tdvnrfK.ijp as ind icated 
in Section ESI Furthermore, iDownes fc Cabri 3 ll2007t) note 
that both these methods can overestimate the flow age and 
thereby underestimate the timescale-based outflow param¬ 
eters. The authors And a more accurate representation of 
tdyn using 1/3 i?iobe/(u), as the intensity-weighted-velocity, 
(v), is probably a better measure of the transverse expansion 
speed of the lobe and where i?iobe is the the perpendicular 
distance from the jet axis (in their jet driven outflows mod¬ 
els). Since the outflows are generally not resolved in this 
perpendicular direction in our data (and in general single 
dish studies), this dynamical-time estimate cannot be tested 
until interferometric observations are obtained. 

Interferometric observations have begun to separate 
som e of the complex regio ns where multiple outflows overlap 
fe.g. iBeuther et ani2002al . Cunningham et al. 2015, submit¬ 
ted to MNRAS), however establishing the source inclina- 
tion to correct tdvn is still difflcult. In li ght of recent work 
dPeters et al.l[2014 . iKlaassen et al.l[2015h even outflows re¬ 
solved on 1-2 arcsecond scales may be a combination of out¬ 
flows from a small cluster of sources (also see Section [5.51) . 
thus single outflow tdyn values may not always be measur¬ 
able. 

The major caveat of the dynamical timescale, on which 
M, P and E depend, is a fairly fundamental one, in that it 
may not a ctually be strongly related to the age of the out¬ 
flows fsee. lCurtis et aLll2010l j. tdyn is model-dependent and 
assumes that material has been accelerated in a shock, prob¬ 
ably by a jet, and is travelling outward at a characteristic 
flow velocity. However, even in the jet scenario, much of the 
high-velocity material is probably turbulently entrained (at 
the bow shock or along the jet sides and outflow cavity walls) 
or accelerated in situ by the passing jet and not have trav¬ 
elled as far as assumed. The outflow may not even trace the 
jet path well at all, and is certainly slower than the jet. Even 
in the case of our example source G078.1224+03.6320 (IRAS 
20126 +4104), the driving je t looks as if it has processed over 
time dShepherd et al.ll2000l j whereas the molecular outflow 
appears to highlight the entirety of its previous and present 



Source Luminosity (Lg) 

Figure 6. The dynamical timescales calculated using the single¬ 
value method, Rmax/i^ma-K, versus the source luminosity. The 
blue- and red-shifted timescales are -^^constant with luminosity, 
there is no significant correlation between the two parameters. 
The open and filled symbols represent YSOs and Hii regions. Only 
sources with D<6kpc are plotted. The dashed an d solid lines are 
the MY SQ and Hll region phase lifetimes from iMottram et alj 

(Eonl). 


path. Calculating fdvn u sing l/3Riobe/(w) as reported by 
iDownes fc Cabrit (l2007l l will alleviate some of these issues 
as this primarily uses outflow parameters at the jet-outflow 
interaction. This dynamical timescale is independent of the 
distance travelled by the outflow from the source, however, 
interferometric observations are a necessity in order to re¬ 
solve the minor axis of the interaction region. 

Additionally, since the jets associated with outflows 
are often observed to leave the natal star-fo rming cores 
and travel f ar into the surrounding diffuse ISM ll Bally et al.l 
I 2 OOI l2012ll . many molecular (CO) outflows may be ‘fos¬ 
sils’ of a past jet event now coasting under momen- 
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turn conservation and their size scale may be set by 
the extent of the core itself. SiO observations could be 
used as an additional tracer to identify more ‘active’ 
outflows, as single dish studies detect board wings for 


outflows, as single disn studies detect board wings tor 
many mid-IR bright sources (e.g. Klaassen fc Wilyir] l2007l : 


iLopez-Sepulcre et al.llioill : iMottram fc Bruntll2012l 'lTRecent 

Sub-Millimeter-Array (SMA) observations by Cunningham 
et al. (2015, submitted to MNRAS) of NCG 2264 indicate 
the strong, collimated SiO (5—4) outflows are only detected 
around the IR-dark, most deeply embedded sources in the 
cluster, which are likely to be the youngest. Such observa¬ 
tions would need to be spatially resolved to ensure the SiO 
emission has the same outflow morp hology of the CO emis¬ 
sion le.g. iDuarte-Cabral et alll2014h . 

The dynamical timescales should therefore not be over¬ 
interpreted as giving accurate source ages or used to calcu¬ 
late accurate accretion rates. They should only be used to es¬ 
timate dynamic parameters for comparison with other stud¬ 
ies following the same methods, and estimating accretion 
rates to an order of magnitude. Our dynamical timescales 
are comparable with the deduced acc retion timescales for 
source s covering a range of masses in iDuarte-Cabral et al.l 
ll2013l) . These authors suggest a roughly constant accretion 
time for the entire mass range, low- to high-mass. If outflows 
are intrinsically related to the accretion phase, this would 
suggest a link between the dynamical timescale of the out¬ 
flows and accretion timescales for the source. Such a rela¬ 
tionship is consistent with what we see in Figure (6] where 
tdyn is ^constant with source luminos i ty (i. e. core mass. 
Paper I). Furthermore, I McKee fc TanI (l2003fl find typical 
timescales for the formation of massive stars as ~ 10® yrs, 
consistent again with our values of tdyn, suggesting tdyn may 
be used as a proxy for accretion timescales. 

That said, it is encouraging to note that the dynami¬ 
cal timescales derived from the outflows compare favourably 
with the MYSO and H ll region phase lifetimes obtained by 
iMottram et al.l ll2011al '). These lifetimes represent the total 
time expected to be spent in each phase, so most sources 
belonging to either category should be younger than these 
values. As shown in Figure [6] the dynamical times for all 
sources are equal to or below the ph ase lifetimes, while also 
older than the yrs predicted bv iDavies et al.l (1201 il l for 
sources to be too faint in the mid-IR to be included in the 
RMS survey. 


5.3 Accretion rates 

In principle, if inclination angles were known and dynamical 
timescales were interpreted as an estimated age of the star 
cluster driving the outflow and an approximate accretion 
time, they could then be used to estimate a time-averaged 
accretion rate (Mace) onto the cluster/core (but, crucially, 
not onto an individual star). 

(Macc) can be crudely estimated using the average tdyn 
for each lobe. In Paper I we showed, from fitting in the 
mass-luminosity plot, that if we assume the luminosity is 
that of an embedded young cluster, then an average ~40 
percent star formation efficiency (SFE) is required in or¬ 
der for the stellar mass of that cluster (with an IMF dis¬ 
tribution of stars) to match the measured core mass. Us¬ 
ing tdyn = R max Iv max ; (Alacc) — Alcore ^ SFEy^tdyn ranges 
from ~1.3xl0“'‘ to ~8.7xl0“® Mq yr“^, with a crude av¬ 


erage of ~2xl0“® Mq yr“^, consiste nt with values used 
in high-mass star formation models (|McKgg__^frajj| 2003; 
Yorke_&_^odcnheimd l2008l : iHosokawa fc Omukail 20011 ; 
Hosokawa et al.l 2010lf . Note, this is a time averaged accre¬ 
tion rate onto a dense cluster core, containing a distribu¬ 
tion of protostars, not onto an individual massive protostar 
alone. 

The luminosity and accretion rate are correlated at 
the P = 0.01 (1 percent, ^2.5 a) level (only using sources 
D<6kpc), although this is driven essentially by the core 
mass itself. Hence the accretion rate and outflow mass are 
also correlated (but at a less significant level P~0.05, or 2 cr; 
see Table O. The main driver being the core mass, which is 
strongly correlated with the outflow mass, suggests that the 
most massive cores are indeed accreting more material and 
therefore the outflow mass itself can be used as a very crude 
proxy (~within an order of magnitude) for the accretion rate 
(under the assumption of tdyn as an order of magnitude age 
estimate). 

Similar arguments, based on the correlation of mass 
flow rates in outflows with source luminosity and, there¬ 
fore, protostellar mass, have been used to support the hy¬ 
pothesis of accelerating accretion rates in the formation of 
massive stars and even to predict birthlines for massive pro¬ 
tostars l|Norberg fc Maederll200fll 'l. Such arguments require 
the assumption that the observed molecular-outflow proper¬ 
ties of MYSOs can be interpreted as relevant to the sequence 
of protostellar evolution, rather than as time-integrated or 
time-averaged quantities produced over the lifetime of the 
forming star, and may therefore be a logical step too far. 


5.4 Mass Flow, Force and Power 

The middle-right and lower panels of Figure [5] show the 
mas s flow rate, force (also ref erred to as momentum flux, 
e.g.. lDuarte-Cabral et al]l2013h and power (commonly called 
mechanical luminosity) of the outflows calculated using the 
spatially varying dynamical timescale. A clear, linear scaling 
is seen in these logarithmic plots, albeit with around an or¬ 
der of magnitude scatter, larger than that seen for the M, P 
and E parameters, due to the range of dynamic timescales at 
a given luminosity (see Section 15.211 . This is reflected in the 
correlation coefficients in Table[5]which are lower than those 
for the latter tdyn-dependent parameters, albeit still highly 
significant. That M, P and E also scale with luminosity is 
usually taken to indicate a common, scalable driving mech¬ 
anism fSection 15.51 discusses the relationship with low-mass 
protostars). 


5.5 Low-Mass Analogues and Clustered Sources 

The scaling of outflow properties between low- and high- 
mass sources could infer a similar driving mechanism, and 
ultimately a similar star formation scenario for high-mass 
stars. We compare our results with outflows from single low- 
mass protostars, known to have jet-driven outflows in their 
early stages of evolution. Our sample of protostellar outflows 
are ema nating from cores asso ciated with massive star for¬ 
mation (iLumsden et al.ll2013l 'l. These cores harbour many 
protostars, some of which are massive, although these cores 
are at varying stages of evolution as identified by the MYSO 
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Table 4. Dynamic timescale and t^yn-dependent parameters calculated for all sources (including those where D>6kpc) with outflows 
or with evidence of outflows, where apertures could be defined. The ‘b’, ‘r’, ‘ave’ and ‘tot’ subscripts indicate the blue-shifted lobe, 
red-shifted lobe, average (of both lobes) and total values of each parameter, fdyn values are In units of 10^ yr, mass flow rates (M) in 
10“^ Mq yr“^, force (P) in 10~® Mq kms“^ yr“^ and power (E) in Lq. The full table is available in the online supplementary material. 
Uncertainties as in Figure [S] are 50 percent for M,P and E, as discussed in the text, uncertainty in how f^yn is calculated and also due 
to inclination correction can make these easily an order of magnitude. 


MSX Source Name 

^dyn,b 

^dyn,r 

^dyn,ave 

Mb 

Mr 

^tot 

Pb 

Pr 

Aot 

Eb 

Er 

Etot 

GOlO.8411-02.5919 

6.0 

7.7 

6.8 

3.8 

4.0 

7.8 

1.2 

1.2 

2.3 

0.4 

0.3 

0.7 

G012.0260-00.0317 

18.6 

23.6 

21.1 

5.3 

3.0 

8.3 

2.8 

1.3 

4.1 

1.6 

0.6 

2.2 

G012.9090-00.2607 

2.5 

1.9 

2.2 

5.2 

28.6 

33.8 

3.4 

26.9 

30.2 

2.4 

26.0 

28.3 

G013.6562-00.5997 

4.3 

7.2 

5.7 

8.3 

1.0 

9.3 

4.0 

0.7 

4.6 

1.9 

0.5 

2.4 

G017.6380+00.1566 

9.8 

8.0 

8.9 

4.7 

8.3 

13.0 

1.8 

4.6 

6.4 

0.7 

2.4 

3.1 

G018.3412+01.7681 

7.7 

6.7 

7.2 

2.6 

3.0 

5.6 

0.5 

1.3 

1.8 

0.1 

0.6 

0.7 

G020.7438-00.0952 

23.2 

27.7 

25.4 

11.9 

3.2 

15.1 

5.9 

1.2 

7.1 

2.7 

0.4 

3.1 

G020.7491-00.0898 

14.7 

15.9 

15.3 

21.3 

6.1 

27.4 

11.7 

2.5 

14.2 

5.8 

0.9 

6.7 

G020.7617-00.0638 

11.3 

4.6 

8.0 

5.7 

13.7 

19.5 

2.6 

7.6 

10.2 

1.1 

4.7 

5.8 

G023.3891+00.1851 

9.9 

8.9 

9.4 

3.4 

4.8 

8.2 

1.3 

2.0 

3.3 

0.6 

1.0 

1.6 


or Hll region classification. It is not necessarily true that 
the massive protostars within these cores are responsible for 
powering the outflows however; i) massive stars may form 
differently to lower mass stars and not power jets that are 
thought to drive outflows; ii) what appears to be a massive 
outflow could be explained by a low-mass protocluster, iii) 
the massive stars in these cores may be too evolved such that 
they no longer power the outflows observed. These scenarios 
are tested below. 


back to the region between Class-0 and Class-I low mass 
outflows (Figure [3 black dotted line). At this stage it is un¬ 
clear whether the shallower slope is due to these cores being 
protoclusters, rather than a sin gle outflow from a single pro¬ 
tostar. As previously noted by BontempsetalJ d 19961) the 
slope of all sources in ICabrit fc Bertoud 1 19921) is shallower 
~ 0.7, although a linear fit (log-log) to only their sources 
where L < 100 Lq indicates a slope of 1. Independent of 
the best fit lines, the outflow force is seen to scale over 
orders of magnitude in luminosity. 


5.5 .1 Scaling of outflow force 

Figure [3 presents the total outflow force of each outflow 
source versus luminosity, for our distance-limited sample 
of massive protostellar cores, together with outflows asso- 
ciated with Clas s -0 and Class-I, young low-mass YSOs from 
Bontemps et al. (Il996l). '^ Class-I, low-mass outflows from 


van der Marel et al.l ( 201^ (their M7 method in determin¬ 


ing force being comparable to ours using P/tdyn, where 
tdyn = Pmax/u) and the pr oposed Class-0 analogue, IR -dark 
high-mass protostars from iDuarte-Cabral et all ll2013l) . For 
consistency, the factor of ~2.9 has been applied to our mas¬ 
sive outflo ws to scale for an aver age 57.3° inclination an¬ 
gle used in iBontemps et ^ (Il996l ). The continuity between 
the low-mass and high-mass samples is striking. When ex- 
tended, the best - fit lin e to the low-mass Class-I sources of 
iBontemps et al.l (1 19961 ). logio F = —5.6 -I- 0.9xlogio ^(Lq), 
intersects directly with our massive outflow sample, and, lies 
slightly below the location of both th e low- and proposed 
high- mass Class-0 analogue sources (iDuarte-Cabral et al.l 
I 2 OI 3 I) . At face value this suggests, at least for lower mass 
protostars, a decrease in outflow force with age. When com¬ 
paring with our more luminous, more massive sources we 
cannot distinguish between MYSO and Hll regions at this 
resolution (Paper I) to establish if they are at a different 
evolutionary stages and there are no clear segregations be¬ 
tween source types observed. Furthermore, offsets between 
different data sets may also mimic such a trend depen¬ 
dent on methodology and inclination corrections. The best 
6 t line to our more massive outflows is slightly shallower, 
Logio P = —4.8 + 0.61x Logio P(L 0 ), and extrapolates 


In terms of outflow force, our IR-bright MYSOs and 
compact Hll regions are positioned as high-mass analogues 
of Class-O/I, low-mass YSOs, if simply extrapolating back 
the fit to our data. A similar correspondence is found if 
we examine the mass-luminosity relationship as in Paper I 
(see Figure El). The high-m ass, Class-0 protostar sample of 
iDuarte-Cabral et al.l (l2013l ) is be positioned below the plot¬ 
ted stellar line as are our cores, albeit with lower luminosi¬ 
ties and masses. Note, the Class O/I classihcation does not 
directly indicate a comparable evolutionary stage (or age) 
for our high-mass sources, as we do not resolve individual 
protostars within the cores. 


The result in Fi gure [3 is very similar to that in 
Ide Villiers et al.l (l2014l ). who examined a sample of distant 
outflows associated with methanol-maser so urces, which are 
reliabl e flags of high-mass star formation llUrauhart et al.l 
l2014bl) . comparing the outflow force to the clump mass. 
This is unsurprising, since their clump masses have a near- 
linear relationship wit h the embedded MYSO luminosity 
dUrauhart et al.l|2014lJ) . as do our cores (Paper I). However, 
given that their sources are much more distant, physically 
larger, more massive and not proven to be a representa¬ 
tive sample, we do not undertake a rigorous comparison. 
Similarly, we acknowledge that there are higher lumi nosity 
targets in the sample accumulated bv lWu et al.l (l2004h with 
which we could also extend our plot, however; some sources 
are much more distant; the sample is not representative; 
sources have rough luminosity estimates; and many differ¬ 
ent methods were used to calculate the outflow parameters. 
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Source Luminosity (Lg) 


Figure 7. Outflow force P versus luminosity for the distance-limited sample (D<6kpc). The filled and open circles represent the force 
calculated with a fixed = RmaxI'Vmax and those summed from the maps where a position variable is used, respectively. 

Both methods p r ovide comparable parameter values. The open and fllled s quares are Class 0 and Clas s I low-mass outflow sources from 
iBontemps et al] lll996h . the filled diamonds are t he '^Class-I sources from Ivan der Marel et alj ll2013lb while the open triangles are the 
high mass outflow Class 0 a nalogue sources f r om IPuarte-Cabral et alJ <l2013l~l . The black solid line is the linear trend of outflow force 
with luminosity presented in lBontemps et alj <ll996l^ for their Class I sources, extended to higher luminosities, while the black dotted 
line is the best fit to our massive outflows extrapolated to lower luminosities. The blue dotted line represent the modelled total outflow 
force from all protostars in a coeval IMF cluster, the blue dashed is the modelled outflow force if only protostars M <8 Mq contribute, 
and the blue dot-dashed line is the modelled outflow force if only the most massive 30 percent of protostars in the cluster provide the 
outflow (see the text for details). 


5.5.2 Are massive protostars driving the outflows? 

When considering the hypothesis of a scaled outflow driv¬ 
ing mechanism, the question arises as to whether there is an 
upper limit to the YSO luminosity range to which this scal¬ 
ing applies, beyond which there might be a different driving 
mechanism and, hence, underlying star-formation process. 
Here we examine whether the most massive stars (M> 8 M©, 
L> 10^ Lq) still form via disc accretion and have jets produc¬ 
ing outflows, or if the massive outflows from luminous cores 
are actually a combination of flows from the many low and 
intermediate-mass sources i.e. can a low/intermediate mass 
protocluster explain the characteristics of a ‘single’ high- 
mass outflow. This also simultaneously tests the scenario 
where the massive protostars in the cluster are more evolved 
and are no longer powering the outflows. Higher-resolution 
observations may eventually provide the answer but, mean¬ 
while, we can use the simple model, as used in Paper I, to 
predict the outflow force from an embedded cluster under 
the assumptions that the protostars have formed coevally 
and that force scales with luminosity. 

We create a population of protostars following the 
Salpeter power-law IMF (using only stellar masses from 0.5 
to 150 Mq) for a range of model cluster masses and then 


calculate the corresponding luminosities of the clu ster us¬ 
ing stellar luminosities from ISalaris fc Cassisil (|200(il . Figure 
5.11) for ma s ses ra nging from 0.5 Mq to 6 Mq and from 
IPavies et al.l (l201ll') for masses >6Mq. The outflow force 
for each protostar in the protocluster is calculated accord- 
ing to the relationship for low-mass, nearby, protostars in 
iBontemps et al.l (llOod l. logio F = —5.6 -I- O.Oxlogio L(Lq), 
which are more likely to be ‘single’ protostars. The blue lines 
in Figure [3 show the results of the luminosities against the 
calculated outflow force in each of our test cases. As noted in 
Paper I, we have not attempted to model a cluster of evolv¬ 
ing protostars in detail, given that high mass stars evolve 
more quickly than lower mass stars, and hence reach their 
ZAMS luminosity more quickly. 

Where the luminosity arises from the complete cluster 
we And the summed total cluster outflow force is closely con¬ 
sistent with the observations (dotted blue line in Figure [7]). 
All protostars up to ~30 Mq must contribute to the outflow 
force in order to match the observations. The M = 30 Mq 
limit is set by the most massive protostar predicted by the 
model in our most luminous core where D< 6 kpc. This result 
supports a similar outflow mechanism for high-mass proto- 
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Table 5. Spearman rank correlation statistics for a range of pa¬ 
rameter relationships. ‘Blue’ and ‘Red’ indicate the blue- and 
red-shifted lobe parameters, respectively; ‘total’ represents the 
total of both lobes except in the case of fdyn 7 correlations, 
where the average value is used. ‘7thick’ ‘7thin’ ^re the Mass- 
spectrum slopes, for which only definite outflow sources with good 
apertures are used in the correlation. All the correlations where 
the P-value is <0.001 are correlated at least at the 0.001 signifi¬ 
cance level. The accretion rate is correlated with source luminos¬ 
ity at the quoted significance level, whereas it is only correlated 
with the outflow mass at the 0.05 level for the quoted p. P-values 
of 0.05, 0.002 and <0.001 represent the ^^2, 3 and >3 (j confidence 
levels. Size is the value indicating the number of sources in each 
correlation. 


Correlation with: 

Blue 


Red 


Total 


(size) 

P 

P-value 

P 

P-value 

P 

P-value 

Lum. with (68): 

Outflow Mass 

0.57 

<0.001 

0.45 

<0.001 

0.55 

<0.001 

Momentum 

0.62 

<0.001 

0.53 

<0.001 

0.61 

<0.001 

Energy 

0.59 

<0.001 

0.54 

<0.001 

0.58 

<0.001 

Mass Flow 

0.54 

<0.001 

0.50 

<0.001 

0.55 

<0.001 

Mech. Force 

0.50 

<0.001 

0.50 

<0.001 

0.51 

<0.001 

Mech. Power. 

0.44 

<0.001 

0.45 

<0.001 

0.45 

<0.001 

Mean Vel. 

0.29 

0.02 

0.34 

0.005 

0.30 

0.12 

Lum. with (32): 

'Ythick 

0.16 

0.39 

0.25 

0.17 

0.31 

0.09 

^thin 

-0.17 

0.35 

-0.27 

0.13 

-0.28 

0.12 

Core Mass 
with (48): 

Outflow Mass 

0.72 

<0.001 

0.64 

<0.001 

0.75 

<0.001 

Momentum 

0.77 

<0.001 

0.74 

<0.001 

0.77 

<0.001 

Energy 

0.70 

<0.001 

0.67 

<0.001 

0.70 

<0.001 

Mass Flow 

0.65 

<0.001 

0.60 

<0.001 

0.62 

<0.001 

Mech. Force 

0.57 

<0.001 

0.52 

<0.001 

0.57 

<0.001 

Mech. Power. 

0.49 

<0.001 

0.45 

0.001 

0.49 

<0.001 

Accretion 

Rate with (48) : 

Outflow Mass 

0.33 

0.02 

0.30 

0.04 

0.37 

0.01 

Source Lum. 





0.43 

0.002 

^dyn with (68): 

Source Lum. 

0.20 

0.11 

0.16 

0.20 

0.17 

0.16 

Momentum 

0.17 

0.15 

0.24 

0.05 

0.19 

0.11 


Outflow 

Energy with (48) : 


Turb. Energy ... ... ... ... 0.73 <0.001 

Binding Energy ... ... ... ... 0.71 <0.001 


stars and thereby supports an upscaled star formation sce¬ 
nario. 

We can also consider the case where massive ‘single’ 
protostars do not form in the same way as their low-mass 
counterparts, and therefore do not produce outflows at all. 
This simultaneously tests whether a low/intermediate proto¬ 
cluster can explain the outflow force, and also if the massive 
protostars have actually stopped powering outflows. We test 
this by assuming only protostars with masses up to ~ 8 M 0 , 
classically low and intermediate mass, contribute to the total 
outflow force. We find a break in the luminosity-force rela¬ 
tion (dashed blue line) at a cluster luminosity of ~6400 Lq. 
This is the protocluster luminosity where the most massive 
protostar reaches Mq . This line can be used to highlight 



Figure 8. Mass-luminosity plot adapted from Paper I. The open 
and filled symbols represent MYSOs and Hll regions respectively, 
except the open tria ngles that are the propose d Class-0 high- 
mass analogues from IPuarte-Cabral et al.l bOl.lh . The two star 
symbols represent Hll regions in the sample that appear to have 
dispersed their core material. The dot-dashed and dashed lines 
represent the luminosity of the most massive star in the IMF clus¬ 
ter and the total luminosity of the cluster re spectively (where the 
star fo rmation efficiency is 50 percent). The lDuarte-Cabral et aU 
boi.-jh Class-0 high-mass analogues are mainly positioned below 
the stellar zero-age-main-sequence (ZAMS) line, which currently 
intersects our cores. This supports the interpretation that our 
cores are high-mass Class-O/I analogues in term of position in a 
mass-luminosity plot. 


two key results. Firstly, the shallower slope of outflow force 
versus luminosity after ~64OOL0 is not consistent with the 
observations. Thus the scenario where massive stars do not 
produce outflows either due to a different formation mech¬ 
anism or being too evolved is unlikely in the case of IMF 
coeval protoclusters with luminosities over ~64OOL0. The 
second result is that, shortwards of ~64OOL0 the line fits 
almost centrally through the observations and at face value 
suggest that low/intermediate mass protoclusters can ex¬ 
plain the outflow force observed in cores where L <6400 L 0 , 
in cases where higher resolution observations are not avail¬ 
able to resolve the sub-structure of the core and identify 
specific outflow drivers. 

The above tests support the interpretation that massive 
protostars power outflows and form in a similar fashion to 
low-mass protostars. As noted above, the most massive pro¬ 
tostars in a cluster are likely more evolved (closer to ZAMS) 
than the low/intermediate mass ones. We attempt to ac¬ 
count for the non-coeval evolution and test the influence 
of the most massive protostars only by curtailing the out¬ 
flow contributing protostars to the most massive 30 percent 
in the protoclusters. The dot-dashed blue line in Figure [7] 
sits slightl y below that ext r apolat ed from the low-mass pro¬ 
tostars in Bontemps et al.l (Il996h but is entirely consistent 


with the data, suggesting that the most massive protostars 
in the clusters are responsible for the outflows alone. Recent 
modelling of outflows from clusters of stars bv IPeters et al.l 
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(|2014) do suggest that the most massive protostars domi¬ 
nate the force and power, even when multiple outflows are 
combined. These authors model four intermediate to mas¬ 
sive protostars that form ronghly coevally (and have simi¬ 
lar masses) from the same accretion disc structure and all 
produce individual outflows with the same common axis of 
projection. On the largest scales, at similar resolution to the 
JCMT observations presented in this paper, a single outflow 
from the system would be seen. Even their sub-arcsecond, 
simulated Atacama Large Millimeter/submillimeter Array 
(ALMA) observations only just begin to separate the bow 
heads of each outflow, but still cannot disentangle the indi¬ 
vidual flows. 


5.5.3 The low-mass scaling relationship 

It is possible that the scalin g relationship assumed above, 
from iBontemps et al.l (Il996ll is not representative, being a 
sample of specifically bipolar low-mass class I sources. In¬ 
deed, other samples appear to fin d shallower slopes e.g. 
F oc (ICabrit fc Bertou3ll992l L However, if we repeat 

onr above analysis with this shallower slope (and placing the 
line snch that it still fits the protostars in the lower-mass 
end, i.e. adjusting the offset), the same key results are ob¬ 
tained. The reported break in luminosity at ~6400 Lq does 
occurs at a lower outflow force, but in general this does not 
change the conclusion that low/intermediate protoclusters 
can be responsible for outflows from cores with luminosities 
below 6400 L 0 . We also reiterate the importance of ensuring 
that different samples of outflow sources use the same analy¬ 
sis methods in order to be fully compatible and thus remove 
the ambiguity of estimating offsets to shift the relationships 
in log-log plots. 

Independent of the relationship adopted from low-mass 
outflows, considering the simplistic model with an ideal, co¬ 
eval Salpeter IMF distribution of protostars, we find the 
most likely conclusion is that vast majority of sources in 
an embedded young cluster contribute to the outflow dy¬ 
namics and that MYSOs with at least ~30 Mq must con¬ 
tribute to the observed flows. The RMS cores are proto¬ 
clusters and most show a single outflow except for a few 
cases where multiple outflows are tentatively seen, with non- 
aligned axes of projection (e.g., GOIO.8411—02.5919 and 
G109.8715-1-02.1156). In these sonrces, the outflow parame¬ 
ters of the system are calculated as a whole (as separation 
of each flow is subjective due to limited spatial resolution) 
and are indistinguishable from the other ‘single’ outflows 
we detect. Hig h -resol ution observations do support this as 
iKlaassen et al. ll2015h present VLA SiO observations where 
three massive protostars transitioning to Hll regions within 
a 2000 an core all have outflows contributing to one large- 
scale (>5000 an) massive outflow. Our model curtailing the 
contributing sources to the most massive 30 percent is likely 
a realistic case for these cores where the massive protostars 
are closer to the ZAMS and dominate over the outflow dy¬ 
namics. These conclusions fully support an upscaled outflow 
mechanism and therefore star formation scenario to massive 
stars. 


Table 6. Linear fit values for luminosity and core mass relation¬ 
ships. The offset and slope are fits in log-log space and as such 
correspond to the equation, Logic (param) = offset + slope X 
LogioL(LQ). The (size) indicates how many sources in the corre¬ 
lation. Only sources with distances <6kpc are used. 


Correlation with: 

(size) 

Offset 

Slope 

Lum. with (68): 

Outflow Mass 

-0.67 ± 0.39 

0.54 ± 0.10 

Momentum 

-0.59 ± 0.42 

0.67 ± 0.10 

Energy 

-0.68 ± 0.52 

0.78 ± 0.13 

Mass Flow 

-5.18 ± 0.34 

0.49 ± 0.08 

Mech. Force 

-5.07 ± 0.46 

0.61 ± 0.11 

Mech. Force inc. corr. 

-4.60 ± 0.46 

0.61 ± 0.11 

Mech. Power. 

-2.92 ± 0.62 

0.72 ± 0.15 

Mean Vel. 

0.41 ± 0.17 

0.12 ± 0.04 


Core Mass 
with (48): 

Outflow Mass 

-0.25 

± 

0.26 

0.77 

± 

0.11 

Momentum 

0.03 

± 

0.25 

0.88 

± 

0.11 

Energy 

0.13 

± 

0.33 

0.98 

± 

0.14 

Mass Flow 

-4.53 

± 

0.23 

0.55 

± 

0.10 

Mech. Force 

-4.24 

± 

0.32 

0.67 

± 

0.13 

Mech. Power. 

-1.93 

± 

0.46 

0.77 

± 

0.20 


5.6 An Alternative Interpretation 

The foregoing discussion explores what might be termed the 
standard model of bipolar molecular outflows, where we in¬ 
terpret the scaling of outflow properties as an indication that 
there is a similar jet driven outflow mechanism for high-mass 
protostars. However, although we see scaling relationships, 
they may not provide any evidence of a similar outflow driv¬ 
ing mechanism, and therefore one cannot assume a scaled up 
star formation scenario. 

An alternative, and considerably simpler explanation 
for the observed correlations in Figure[5]for these high-mass 
outflows yet to be considered, is that the fnndamental rela¬ 
tionship is that between core m ass and outflow mass, a s seen 
in Figure[9]and also reported in ide Villiers et ^ I (l2014h . and 
that all the relationships between the tuyn-dependent out¬ 
flow properties and the source luminosity stem from this. 
There are several pieces of evidence in the data that might 
point to this model. Firstly, the correlation coefficient for 
the relationship between outflow mass and core mass, along 
with that of outflow momentum versus core mass (see Table 
O, is the highest of those measured and therefore contains 
the least scatter, suggesting that one of these two may con¬ 
stitute the basic relationship. 

Of the other calculated dynamical outflow parameters, 
P and E depend on v and v^, while P and E are calculated 
from Mv/tdyn ~ Mv^/R and 0.5Mu^/tdyn ~ Mv^/R and 
so, in addition to the accelerated mass in the outflow, de¬ 
pend only on various powers of the flow velocity and on the 
outflow lobe size. Figure | 6 ] shows that tdyn = Rmaxivmax is 
independent of luminosity, while the evidence for a correla¬ 
tion between the weighted-mean velocity and luminosity is 
only barely significant in this sample (at the 2—3a level, see 
Figure [To] and Table[5ll. Note that, since the momentum and 
energy are calculated via the summation over every veloc¬ 
ity bin and spatial position, the velocity used is the ‘bulk’ 
intensity-weighted velocity (v). The slope of a least-squares 
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Figure 9. Core mass plotted against the optical-depth-corrected 
outflow mass for sources whose core masses are flagged as good 
(see Paper I). There is a clear correlation between the two quanti¬ 
ties. The open and filled symbols represent YSOs and Hii regions 
(as Paper I sources that are classified as YSO/Hii are identifled as 
Hii here), the circles and squares are for ‘Y’ and ‘S’ aperture flags, 
respectively. The sources circled are those at a distance greater 
than 6 kpc while those with a cross symbol are sources with some 
outflow evidence (flagged ‘M’). A 50-percent error is plotted for 
the core mass determination whereas that for outflow masses is 
the average calculated error as Figure [5] 



Source Luminosity (Lg) 


Figure 10. Plot of the average intensity-weighted-mean, ‘bulk’, 
outflow velocity, calculated via {Px,ylMx,y)-, against source lumi¬ 
nosity. There is only a weak correlation over the large range of 
luminosity, although note each outflow has a range of velocities as 
used in the calculation of P and E. The open and filled symbols 
represent YSOs and Hii regions, whereas the circles and squares 
are for ‘Y’ and ‘S’ aperture flags respectively. 
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linear fit (log-log) to the (v) vs L relationship is ^0.12±0.04, 
see Tabled 

With this in mind, the second piece of evidence in the 
data is that the details of the relationships with lumiirosity 
seen in Figure [5] are consistent with the way in which the 
dynamical time based properties are calculated. The slope of 
the linear fit (log-log) to the dependence of the outflow mass 
on L is 0.54±0.10, that of the flow momentum is 0.67±0.10 
and of flow energy 0.78±0.13 and the relatioirships involviirg 
the mass flow rate, force and power are similar. These slopes 
differ by only Icr but are consistent with the sequence of 
increasing powers of v in their calculation (refer to Table 
[6]). Importantly, we see that the scatter in the relationships 
also increases in the same sequeirce. This suggests that the 
additional powers of v iir calculatiirg P and E, etc., is not 
revealing a more fundamental physical relationship and that 
A/fiow has the most basic relationship with L, although this, 
in turn arises because L is linearly dependent on Mcore. 

The hypothesis is, therefore, that the entrained mass is 
the fundamental property of the outflow, that this depends 
only on the mass available in the dense core, and that the 
core mass is the only parameter amongst those examined 
here with a direct physical relationship to t he luminosity of 
the em bedded protostars (see Paper I and lUrauhart et al.l 
l2014bl '). Thus, all the correlations between the flow, force 
and power with luminosity could arise simply because there 
is more mass available to be entrained in massive cores, and 
these massive cores themselves also tend to contain more 
luminous protostars, i.e. ranging from single sources, to a 
cluster, the entrained mass of the outflow is directly related 
to the mass available to each individual protostar. Such that 
in the case of a cluster, each protostar drives an outflow re¬ 
lated to the protostar’s envelope mass, being a fraction of 
the total core mass fe.g. iBontemps et al.lll99B b Such an in¬ 
terpretation, of course, renders any conclusion regarding the 
outflow-driving mechanism suspect, since this remains un¬ 
constrained and any process that accelerates the core gas 
with similar entrainment efficiencies will produce essentially 
the same effect. This, in turn, implies that accretion rates 
may not be directly inf erable from mass o u tflow rates for 
high-m a ss sources as per Richer et al.l ll2000ll ; iBeuther et al.l 


(l2002df : Ide Villiers et al 


ll2014l'l as these assume a partic¬ 


ular family of momentum-conserving driving mechanisms 
with particular efficiencies. Indeed, even if more massive 
protostars drive more powerful outflows, and have higher 
infall rates (which are intrinsically linked to accretion), one 
would require knowledge of the efficiencies of the outflow 
mechanism for high-mass protostars before inferring accre¬ 
tion rates. 

There are several a priori arguments that may support 
this hypothesis. Firstly, r egard less of how the outflow mate¬ 
rial is entrained, |d^ 3 (Il984ll pointed out that the mech¬ 
anism accelerating the gas of the molecular outflow should 
be either momentum- or energy-conserving, but not both. 
Therefore, the virtually equal correlations of outflow force 
and outflow power with source luminosity always found in 
studies such as this cannot simultaneously represent a phys¬ 
ical relationship. The terminal wind/jet velocity is required 
as an input for calculation of the outflow force or outf low lu¬ 
mino sity in the different regimes (see eqns. 25 and 26. lDvsonl 
Il984h . Using the wind/jet speed should indicate a constant 
value of either outflow force or outflow luminosity, while the 
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other parameter scales with source luminosity, dependent on 
the nature of the outflow (i.e. either momentum- or energy- 
conserving) . 

Secondly, in low-mass sources, the jets are known 
to leave the densest r egions of the molecular cloud (e.g. 
iBallv et al.ll200|2012| ~) and travel much further out into the 
diffuse interstellar environment (often many parsecs but, in 
any case, on a scale much larger than the molecular out¬ 
flow), only showing up as bow-shocks when stopped by the 
diffuse ISM. Thus, an unknown fraction of the jet energy or 
momentum, depending on what is conserved — most likely 
momentum — must be deposited in the molecular outflow 
component and this should create enough scatter in the mea¬ 
sured parameters of the latter to wipe out any correlations. 
Yet low-mass out flow sources produce t he most significant 
correlations (e.g. iBontemos et al.l Il99fil '). Whether or not 
the outflows from high-mass sources are mostly jet-driven 
is still uncertain (e.g. as the case for GOIO.8411—02.5919, 
G078.1224-b03.6320 and G109.8715-b02.1156, see above). 

Furthermore, the parameters of the la rge-scale Orion 
molecular outflow (e.g. lErickson et al.l[l982l i are consistent 
with t he standard dynamical relationships dBallv fc Ladal 
1 19831 : ICabrit fc Bertou3 Il992l 'l but this outflow is now 
thought to have been caused by a single impulsive, po ssibly 
explosive event l|ZaData et al.l I 2 OO 9 I : iBallv et al.l 1201 il l , un¬ 
like the jet-driven flows seen in at least the low-luminosity 
sources. This alone must raise doubts over the conclusion 
that the universal correlations imply a single outflow-driving 
mechanism. 

Of course, this mass-only model also predicts that the 
outflow length R will be determined mainly by the density 
distribution and size of the core. This does not appear to 
be the case (Figure fTTI) . however it may be that such a re¬ 
lationship is masked by scatter from the random outflow 
orientation and by highly variable collimation. 

Another consequence of this model is that the dynami¬ 
cal timescales obtained from R/v cannot be related to real 
source ages, except in the youngest flows in which a driving 
jet has not yet reached the edge of the core. In sources in 
which a jet has left the dense core, tdyn may give only an 
estimate of the time taken for the jet driven shock to reach 
the low-density outer regions of the core envelope. 

The power law relationship in Figure [9] is Moutflow oc 
This is not significantly different enough from a 
direct linear relationship to warrant too much explanation. 
If it were, then there would have to be some relationship be¬ 
tween the core mass and the fraction of the core affected by 
the outflow (diminishing if flatter than linear, e.g.) implying 
some kind of mass-dependent collimation of the flow. 

If the preceding model survives further analysis, then 
the collimated nature of massive protostellar outflows, i.e., 
their bipolarity, is the remaining feature that unites them 
to the low-mass version of the phenomenon. This may be 
due to the ubiquitous presence of disc structures related to 
accretion. Bipolarity is clear for some of the outflows from 
this single-dish survey. Some synthetic observations repro¬ 
duce the bipolar morphology of high-m ass outflows, even if 
driven by ionisation feedback alone (e.g. IPeters et al1l2012ll . 
However, the predicted outflow parameter values are lower 
than those observed and hence magnetic forces are again 
suggested as an acceleratio n mechanism. 

The recent review by iLi et al.l (l2014l ~) details that, ir- 



Core Radius (pc) 

Figure 11. Umax for blue- and red-shifted outflow lobes ver¬ 
sus the core radius from Paper I. Open and filled symbols are 
MYSOs and Hll regions, respectively, while the colours represent 
the blue- and red-shifted lobes, and circles and squares are ‘Y’ 
and ‘S’ aperture flags. Clearly there is no correlation between 
these parameters. 


respective of the numerical modelling approach, successful 
models of jets all include some aspect of a rotating disc and a 
magnetic collimation and driving force. Only a rare few high- 
mass sources have been independently studied and indicate a 
sc aled-up picture of l ow-mass star format ion (e.g. Cep A HW 
2; IPatel et al.l l2005l . rotatin g material, |Curie l_eL_^ I 2 OO 6 I . 
proper motion radio jets, and IVlemmings et al. II 2 OIOI outflow 
and jet align e d ma gnetic fields). As argued in the review by 
[Richer et al.l (l200Cllj a high-resolution, high-sensitivity study 
of a significant sample of massive protostellar outflows is re¬ 
quired. Such a study is still lacking, although the sample in 
this work provides a firm footing for such an investigation as 
these high-mass star forming regions clearly drive outflows. 


5.7 The Mass-Velocity Relationship 

In the earlier stages of massive star formation, the convective 
surface layers of the protostar may be conducive to the pro¬ 
duction of magnetic fields which then supply the magneto¬ 
hydrodynamic forces r equired for jet produc tion, powered 
by the stellar rotation dHosokawa et al.ll20ldi . Jets, driven 
and regulated via the complex relationship of MHD forces 
and disc accreti on, may theref ore be inherent to massive-star 
formation (see, iLi et al.ll2014l i. 

One method used to investigate jets as the powering 
mechanism of molecular outflows is to examine the mass 
spectrum, i.e., the power-law relation s hip between inten¬ 
sity and velocity llLada fc Fichl 1 19961 : [Richer et al.l [2000[ : 
[Ridge fc Moore|[200lh . in which a break can often be identi¬ 
fied in the higher-velocity flow that may be due to molecular 
disso ci ation caused by a jet s hock ([Downes fc Cabrig [20031 . 
[ 2 OO 7 II . [Ridge fc Moore! ([2001] ') note that optical-depth cor¬ 
rection must be undertaken to establish an accurate slope for 
the relationship, especially at lower velocities. The left and 
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Figure 12. Left and Right: Mass-spectrum slopes for all sources plotted against luminosity for the optical-depth-corrected velocities 
and the optically thin, higher-velocity ranges in the emission-line wings. The slopes of both velocity regions are comparable. Open and 
filled symbols are MYSOs and Hll regions, respectively, while the colours represent the blue- and red-shifted lobes. 


right panels of Figure [12] show the results for both optical- 
depth-corrected and optically thin velocity regimes. 

Generally, the sl opes are consistent with those found 
for lo w-mass outflows l|Richer et al.ll200(il : IStoiimirovic et al.l 
l2006l 'l. However, no clear breaks are seen in the higher- 
velocity, optically thin regions and single linear relation¬ 
ships fit the mass-spectra of the individual cores reason¬ 
ably well. The optical-depth-corrected and optically thin 
slopes also span the same ranges. There is no evidence 
of a recently accel er ated or coasting out flow component 
dRicher et al.l [20oOi . IPlunkett et al.l I 2 OI 5 I ') do not report 
breaks in the slope of the combined outflows from the two 
clusters they investigate. This confirms that we could not 
expect to see a break in any of our high-mass cores as they 
likely contain multiple, possibly different direction and in¬ 
clination outflows, which when combined would wash-out 
any underlying breaks that are seen for single, jet-driven 
low-mass outflows. The similarity of the slopes for low- 
and high-mass sources could in itself, be interpreted as due 
to a similar outflow-acceleration mechanism. Higher sen¬ 
sitivity and higher resolution observations are required to 
detect the very high velocity flow and con f irm i ts origin 
from a single driving source. iLebron et al.l (l2006l 'l note a 
steepening mass-spectrum for the high velocity outflow > 
40 kms“^ in their high-sensitivity single dish observations of 
G078.1224-1-03.6320, whereas the JCMT observations here 
only detect a Umax of ~40kms“^. 


5.8 Impact on the Natal Core 


Simulations and observations suggest that molecular out¬ 
flows both can, and cannot have a profound effect on the 
natal core (e.g. Arce et al. ] Loid: iM ottram fc Brunj I 2 OI 2 I : 


iFederrath et al.ll2014l and also see lLi et al.ll2014l 'l. data 

presented in Paper I show 30 sources (D<6 kpc) have smooth 
velocity gradients across the cores, six of which are aligned 


with the outflows presented here, suggestive of their strong 
influence. We can compare the energy content of the outflows 
to the turbulent kinetic energy in the core gas and the grav- 
itational binding ene rgy of the core (assuming p{r) oc r~^, 
IShenherd et al.l2007l 'l. In the case where the thermal motions 
contribute little to the measured full width half maximum 
(FWHM) of the core line emission, the turbulent energy 
can be calculated via Eturb = (3/16 In 2)Mcore xFWHM^ 
dArce fc Goodm^l200ll l. Figure [13] shows the turbulent en¬ 
ergy, calculated using Mcore from Paper I and where the 
velocity FWHM is measured using C^®0 (3—2), versus the 
outflow energy. For the binding energy calculation, the core 
radii listed in Paper I are used. 


The outflow energy on average can fully account for the 
turbulent energy of the core and also equates to ~85 percent 
of the binding energy, even without inclination-angle correc¬ 
tion. Hence we might conclude that outflows from massive 
YSOs contribute significantly to the core turbulent kinetic 
energy and are a significa nt source of mechanical feedback. 
I Cunningham et al.l ( 200d l conclude that jet-driven outflows 
act to maintain the turbulence in a molecular cloud, pro¬ 
vided there was already an initial disruption. However, Fig¬ 
ure M (right), shows that the cores that do not contain 
outflows have turbulent kinetic energies that are consistent 
with those that do. This suggests that the core turbulence 
(on scales of ~ 0.4 to ~ 2 pc) is not driven by the local input 
from outflows. If the jets that accelerate the molecular flows 
tend to leave the dense cores, only a fraction of their mo¬ 
mentum and energy might be deposited into the dense core 
gas, and this only in local re gions. In the low-mass core, B59, 
iDuarte-Cabral et al.l (I 2 OI 2 I ') see remnant U-shaped cavities 
and ridges that are a result of the direct impact of outflows 
on the less dense, local material. However, they also see other 
velocity structures in the C^®0 data of the denser core, e.g. 
gradients and infall motions, and although the outflows do 
have enough energy to fully drive the turbulence, it is not 
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Figure 13. Left: Turbulent energy of the cores plotted against the outflow energy. The open circles represent sources where a core 
mass is established while the filled ones are subset where core masses are flag 0, 1, and 2 in Paper I, i.e. reliable estimates. On average 
without inclination correction the outflow energy can fully account for the turbulent energy. Right: Turbulent energy plotted against 
source luminosity. The open and closed circles are the same as the left figure, the squares are for cores without outflows (open and closed 
have the same meaning as the circles). It is clear that cores without outflows appear to have the same range of turbulent energy values. 


straightforward to conclude that the outflows alone are the 
predominant production mechanism. 

For our cores there are many plausible explications how¬ 
ever, simply the outflows are not detected in these sources 
(e.g. plane of sky confusion, or weak) or there are al¬ 
ternative inputs of turbulence (local Hll regions, external 
winds/shocks, or interactions of the sources within the cores 
below our resolution). Since turbulent energy tends to flow 
from larger to smaller scales, we do not expect such localised 
phenomena to transfer energy to the whole cloud for exam¬ 
ple. Although, if energy is transferred on core scales at a 
cavity wall interface, we would require higher spatial resolu¬ 
tion to isolate these i nteraction regions where th e cavity is 
within the beam (e.g. iDuarte-Cabral et al.ll2012) . 


6 SUMMARY 

From a sample of 99 sources, 65 and 20 have been identified 
to be definitely driving outflows and have some evidence for 
outflow, respectively. The remaining 14 show no signs of out¬ 
flows and have Gaussian line profiles. For a distance- 

limited subset (D<6kpc), 59 have definite outflows, 17 have 
some evidence of outflows and 13 do not have evidence of 
outflows. 

The kinematic and dynamic parameters have been cal¬ 
culated for all sources with outflows and with evidence for 
outflows, but only where it was clear how to separate the 
outflow material from the diffuse ambient emission in inte¬ 
grated maps. Furthermore, M, P and E calculated using 
a position-variable dynamical timescale tdyn(a;,!/) (applied to 
massive protostellar outflows for the first time) are shown 
to be consistent with those obtained using a single tdyn es¬ 
timate via Rmax/umax. Using these dynamical timescales 


and SFE values established in Paper I, we infer a time- 
averaged accretion rate of ^2x10 “® Mq onto these massive 
star-forming cores. 

All outflow parameters scale with source luminosity 
and core mass. There are no intrinsic differences between 
the cores classified as MYSOs or Hll regions. The out¬ 
flows are driven by many protostars within a cluster (given 
the resolution of the observations). Specifically, the out¬ 
flow force scales directly from samples of low-mass proto¬ 
stars, suggesting that our sample consists of high-mass Class 
O/I analogues, in terms of position in luminosity-outflow 
force, and luminosity-mass diagrams. The relationship be¬ 
tween outflow force and luminosity is consistent with all 
the sources in a coeval star-forming cluster contributing to 
the outflow force. Models in which the massive protostars 
do not contribute to the outflows are not consistent with 
the observations, although these models do suggest outflows 
from low/intermediate mass protoclusters can explain ‘sin¬ 
gle’ outflows from cores with luminosities <6400 Lq. Models 
curtailing the mass of the contributing protostars to only the 
most massive 30 percent are coincident with the observa¬ 
tions, suggesting that the massive protostars dominate over 
the outflow dynamics. The data support a scaling up of the 
star formation process for massive protostars up to ~30 Mq 
in a star-forming cluster. 

An alternative interpretation of the scaling relationships 
is that the molecular outflow parameters are determined al¬ 
most entirely by the entrained mass, which is set by the 
available core mass. This model does not require a single 
driving mechanism for the molecular outflow as any form 
of mass acceleration, explosive or continuous, would pro¬ 
duce the observed correlations. The bipolar nature of some 
of these high-mass outflows may be the key to supporting a 
scalable driving mechanism from low-mass sources. 
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Although the outflow energetics are comparable to the 
turbulent energy iir the deirse cores, we And that cores with 
outflows have similar kinetic energies to those without. Such 
turbulence could be provided by alternative means or by un¬ 
detected outflows. At the scales probed we cannot establish 
the local impact these outflows have. 

The sample presented is ideal for high-resolution, high 
sensitivity follow-ups in order to disentangle the smaller 
groups of protostars driving the outflows, establish their 
bipolar nature and the underlying driving mechanism. 
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APPENDIX A: COLUMN DENSITY AND MASS CALCULATION 

In this appendix the column density and mass equations are derived following from the result of I Garden et alj lll99lli . except 
for the CO(3—2) transition. The total column density of a linear, rigid rotor molecule under conditions of local thermodynamic 
equilibrium (LTE), with the populations of all levels characterised by a single excitation temperature, Tex, is obtained from 
the integral of the optical depth over the line profile: 


_ 3k exp[hBJ{J + l)/fcrex] Tex + hB/3k 
8-K^BfB (J + 1) [1 — exp(—hi//fcTex)] 

where B is the rotational constant, /r is the permanent dipole moment of the molecule and J is the rotational quantum 
number of the lower state, in this case J = 2 for the CO(3—2) transition, k and h are the Boltzmann and Planck constants 
respectively. The excitation temperature Tf.^ is solved for, as s hown in Equation [21 Section [3.II 
Here the approximation for f Tv dv follows feuckle et al.l ll2niflh for the case where r 




, dv 
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Tu dv = 


hv / 1 

k y exp(/iiz/fcrex) — 1 


1 

exp(/l!z/fcTcmb) - 1 


-1 
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[1 - exp(-T)] 


j Tmb dv 


(A2) 


The brightness temperature, Tmb, is the antenna temperature of the telescope divided by the beam efficiency, TX/r^mb, 
and corresponds to the Rayleigh-Jeans brightness of a source minus the brightness of the cosmic microwave background with 
temperature, Tcmb = 2.73 K, over the beam. Combining equations im and IA2I in the limit where Tex ^ Tcmb results in the 
column density: 


_ 3k exp[/iBJ(J + l)/fcTex] 1 Te^ + hB/3k f 1 

^ (J+1) (hiz/fc) [exp(-ht//fcTex)] y [1 - exp(-r)] ^ ’ 

where the permanent dipole moment for *^^00 is 0.1098 Debye for llChackerian fc Tippind [198311 . Conforming to cgs 

units typically used in such analysis i3=58.14GHz, fc=1.381x 10~^® ergK~^, h=6.626x 10~^^ ergs, t/(^^CO) = 345.79599 GHz, 
velocity v is in kms“^, /r(^^CO) = 0.1098 x 10“^® StatCcm (where IstatC = 1 g^^^ cm®^^ s^ = 1 erg^^® cm^^®), r becomes 
ri 2 , the calculated optical depth of the ®^CO line ('Section l3.ll) and Tex is the calculated excitation temperature. The column 
density for the (3—2) transition is therefore: 


A(^®CO) = 4.78 X 10 


12 exp(16.74/Tex) (Tex + 0.93) 
exp(-16.59/Tex) 


/ 


Tmb 


T12 


[1 - exp(-ri 2 )] 


dv , cm 


The mass can then be calculated directly from the column density via: 


(A4) 


(A5) 

where /ig = 1.36 is the total gas mass relative to H 2 , the abundance ratio H 2 /^^CO = lO'^, and D is the distance of the 
source to the Sun, in kpc. D is the solid angle corresponding to the emission in one pixel of the maps used in this work. Thus 
including the conversion factors the core gas mass in solar masses (Mq) is calculated for every pixel of the outflow lobe maps 
using Equation IA6I The total outflow lobe masses as reported in Table [2] is the summation of the pixel masses within the 
defined outflow apertures. 
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